library(here)
## here() starts at /home/n.brouwer/MEP
library(ggplot2)
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2
## ──
## ✔ tibble  3.1.8     ✔ dplyr   1.1.0
## ✔ tidyr   1.3.0     ✔ stringr 1.5.0
## ✔ readr   2.1.3     ✔ forcats 1.0.0
## ✔ purrr   1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(ggrepel)
library(fst)
source(here("UtilityFunctions.R"))
library(Hmisc)
## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
## 
## Attaching package: 'Hmisc'
## 
## The following objects are masked from 'package:dplyr':
## 
##     src, summarize
## 
## The following objects are masked from 'package:base':
## 
##     format.pval, units
library(knitr)
library(tidyverse)
library(corrplot)
## corrplot 0.92 loaded
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine

Load data

combinations <- readRDS(here('scratch/combinations_selection.rds'))

get_parameters <- function(combinations){
  params <- c()
  for (c in combinations){
    tryCatch({
    success_model <- readRDS(here(paste('scratch/success_models/success_models_', c,'.rds', sep='')))
    params <- append(params, list(success_model[[3]]))
    names(params)[length(params)] <- c
    },error=function(e){print(e)}, warning=function(w){print(w)})
    
  }
  
  params_result <- bind_rows(params, .id='phenotype_combo')
  return(params_result)
}

# Remove outliers
all_parameters <- get_parameters(combinations) %>%
  separate(phenotype_combo, into=c('phenotype_from','phenotype_to'),sep='_to_', remove = FALSE) %>%
  filter(a > 0.5) %>%
  filter(b > 8) %>%
  rename(shape =a , scale = b)
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_B cells_to_CD57^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_Ki67^{+}_to_CK^{lo}ER^{med}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_MHC^{hi}CD15^{+}_to_B cells.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_Fibroblasts FSP1^{+}_to_MHC^{hi}CD15^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_MHC^{hi}CD15^{+}_to_Fibroblasts FSP1^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_Ep CD57^{+}_to_Macrophages & granulocytes.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_Macrophages & granulocytes_to_Ep CD57^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_MHC I^{hi}CD57^{+}_to_MHC^{hi}CD15^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_MHC^{hi}CD15^{+}_to_MHC I^{hi}CD57^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_Fibroblasts FSP1^{+}_to_Macrophages & granulocytes.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_CK8-18^{+} ER^{hi}_to_HER2^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_Fibroblasts FSP1^{+}_to_CD4^{+} T cells.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_Fibroblasts FSP1^{+}_to_HER2^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_CD4^{+} T cells & APCs_to_CD57^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_CD57^{+}_to_CD4^{+} T cells & APCs.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_MHC I^{hi}CD57^{+}_to_Ep Ki67^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_CD57^{+}_to_ER^{hi}CXCL12^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_CD57^{+}_to_Granulocytes.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_Granulocytes_to_CD57^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_HER2^{+}_to_MHC I^{hi}CD57^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_MHC I^{hi}CD57^{+}_to_HER2^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_Ep CD57^{+}_to_MHC^{hi}CD15^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_MHC^{hi}CD15^{+}_to_Ep CD57^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_CK^{lo}ER^{med}_to_Granulocytes.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_Granulocytes_to_CK^{lo}ER^{med}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_Granulocytes_to_MHC I^{hi}CD57^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_MHC I^{hi}CD57^{+}_to_Granulocytes.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_CD57^{+}_to_HER2^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_HER2^{+}_to_CD57^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_CK8-18^{hi}ER^{lo}_to_HER2^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_HER2^{+}_to_CK8-18^{hi}ER^{lo}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_CK8-18^{hi}ER^{lo}_to_Ki67^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_Ki67^{+}_to_CK8-18^{hi}ER^{lo}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_T_{Reg} & T_{Ex}_to_CK8-18^{hi}ER^{lo}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_CK^{+} CXCL12^{+}_to_MHC^{hi}CD15^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_MHC^{hi}CD15^{+}_to_CK^{+} CXCL12^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_CK8-18^{+} ER^{hi}_to_MHC^{hi}CD15^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_MHC^{hi}CD15^{+}_to_CK8-18^{+} ER^{hi}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_ER^{hi}CXCL12^{+}_to_MHC^{hi}CD15^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_MHC^{hi}CD15^{+}_to_ER^{hi}CXCL12^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_CK8-18^{hi}ER^{lo}_to_MHC^{hi}CD15^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_MHC^{hi}CD15^{+}_to_CK8-18^{hi}ER^{lo}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_CD57^{+}_to_MHC^{hi}CD15^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_MHC^{hi}CD15^{+}_to_CD57^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_CD57^{+}_to_CK8-18^{+} ER^{hi}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_CK8-18^{+} ER^{hi}_to_CD57^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_CK8-18^{+} ER^{hi}_to_Granulocytes.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_Granulocytes_to_CK8-18^{+} ER^{hi}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_CD57^{+}_to_Fibroblasts FSP1^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_Fibroblasts FSP1^{+}_to_CD57^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_CD57^{+}_to_Ep Ki67^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_CD57^{+}_to_CK8-18^{hi}CXCL12^{hi}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_CK8-18^{hi}CXCL12^{hi}_to_CD57^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_Fibroblasts FSP1^{+}_to_MHC I^{hi}CD57^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_MHC I^{hi}CD57^{+}_to_Fibroblasts FSP1^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_CD57^{+}_to_CK^{+} CXCL12^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_Fibroblasts FSP1^{+}_to_MHC I & II^{hi}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_Ki67^{+}_to_Basal.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_Ep CD57^{+}_to_Ki67^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_Ki67^{+}_to_Ep CD57^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_CD57^{+}_to_CK8-18^{hi}ER^{lo}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_CK8-18^{hi}ER^{lo}_to_CD57^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_Macrophages_to_CD57^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_CK8-18^{hi}ER^{lo}_to_MHC I^{hi}CD57^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_MHC I^{hi}CD57^{+}_to_CK8-18^{hi}ER^{lo}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_CD57^{+}_to_Macrophages & granulocytes.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_Macrophages & granulocytes_to_CD57^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_CD57^{+}_to_Ki67^{+}.rds', probable reason 'No such file or directory'>
## <simpleWarning in gzfile(file, "rb"): cannot open compressed file '/home/n.brouwer/MEP/scratch/success_models/success_models_Ki67^{+}_to_CD57^{+}.rds', probable reason 'No such file or directory'>
not_estimated_combinations <- setdiff(combinations, unique(all_parameters$phenotype_combo))

# Compute means of each combination
means <- all_parameters %>% 
  group_by(phenotype_from, phenotype_to ) %>% 
  summarise(mean_shape = mean(shape),
            mean_scale  = mean(scale))
## `summarise()` has grouped output by 'phenotype_from'. You can override using
## the `.groups` argument.
all_distances_data <- readRDS(file  = here('scratch/AUCScaledSlides_300_ALL.rds'))
all_distances_data <- all_distances_data %>% as_tibble %>%
  mutate(distance_window = WinMean) %>%
  mutate(phenotype_combo = paste(phenotype_from, phenotype_to, sep='_to_')) %>%
  dplyr::select(tnumber, phenotype_combo, `N.per.mm2.scaled`, distance_window) %>%
  mutate(new = `N.per.mm2.scaled` * 1000) %>%
  mutate(new =round(new))

Failed parameter estimations

Which combinations could not be estimated?

initial_parameters <- tibble()
for (c in not_estimated_combinations){
    init_params <- readRDS(here(paste('scratch/init_parameters/init_params_', c,'.rds', sep='')))
    initial_parameters <- bind_rows(initial_parameters, init_params)
}

initial_parameters <- merge(x= initial_parameters %>%
                              filter(term=='shape') %>%
                              select(c(tnumber,combo, estimate))%>%
                              rename(shape = estimate),
                            y= initial_parameters %>%
                              filter(term=='scale') %>%
                              select(c(tnumber,combo, estimate)) %>% 
                              rename(scale = estimate),
                            by.x=c('tnumber', 'combo'), 
                            by.y= c('tnumber', 'combo')) %>% 
                      separate(combo, into=c('phenotype_from','phenotype_to'),sep='_to_', remove=FALSE)

print(str_sort(not_estimated_combinations))
##  [1] "B cells_to_CD57^{+}"                               
##  [2] "CD4^{+} T cells & APCs_to_CD57^{+}"                
##  [3] "CD57^{+}_to_CD4^{+} T cells & APCs"                
##  [4] "CD57^{+}_to_CK^{+} CXCL12^{+}"                     
##  [5] "CD57^{+}_to_CK8-18^{+} ER^{hi}"                    
##  [6] "CD57^{+}_to_CK8-18^{hi}CXCL12^{hi}"                
##  [7] "CD57^{+}_to_CK8-18^{hi}ER^{lo}"                    
##  [8] "CD57^{+}_to_Ep Ki67^{+}"                           
##  [9] "CD57^{+}_to_ER^{hi}CXCL12^{+}"                     
## [10] "CD57^{+}_to_Fibroblasts FSP1^{+}"                  
## [11] "CD57^{+}_to_Granulocytes"                          
## [12] "CD57^{+}_to_HER2^{+}"                              
## [13] "CD57^{+}_to_Ki67^{+}"                              
## [14] "CD57^{+}_to_Macrophages & granulocytes"            
## [15] "CD57^{+}_to_MHC^{hi}CD15^{+}"                      
## [16] "CK^{+} CXCL12^{+}_to_MHC^{hi}CD15^{+}"             
## [17] "CK^{lo}ER^{med}_to_Granulocytes"                   
## [18] "CK8-18^{+} ER^{hi}_to_CD57^{+}"                    
## [19] "CK8-18^{+} ER^{hi}_to_Granulocytes"                
## [20] "CK8-18^{+} ER^{hi}_to_HER2^{+}"                    
## [21] "CK8-18^{+} ER^{hi}_to_MHC^{hi}CD15^{+}"            
## [22] "CK8-18^{hi}CXCL12^{hi}_to_CD57^{+}"                
## [23] "CK8-18^{hi}ER^{lo}_to_CD57^{+}"                    
## [24] "CK8-18^{hi}ER^{lo}_to_HER2^{+}"                    
## [25] "CK8-18^{hi}ER^{lo}_to_Ki67^{+}"                    
## [26] "CK8-18^{hi}ER^{lo}_to_MHC I^{hi}CD57^{+}"          
## [27] "CK8-18^{hi}ER^{lo}_to_MHC^{hi}CD15^{+}"            
## [28] "Ep CD57^{+}_to_Ki67^{+}"                           
## [29] "Ep CD57^{+}_to_Macrophages & granulocytes"         
## [30] "Ep CD57^{+}_to_MHC^{hi}CD15^{+}"                   
## [31] "ER^{hi}CXCL12^{+}_to_MHC^{hi}CD15^{+}"             
## [32] "Fibroblasts FSP1^{+}_to_CD4^{+} T cells"           
## [33] "Fibroblasts FSP1^{+}_to_CD57^{+}"                  
## [34] "Fibroblasts FSP1^{+}_to_HER2^{+}"                  
## [35] "Fibroblasts FSP1^{+}_to_Macrophages & granulocytes"
## [36] "Fibroblasts FSP1^{+}_to_MHC I & II^{hi}"           
## [37] "Fibroblasts FSP1^{+}_to_MHC I^{hi}CD57^{+}"        
## [38] "Fibroblasts FSP1^{+}_to_MHC^{hi}CD15^{+}"          
## [39] "Granulocytes_to_CD57^{+}"                          
## [40] "Granulocytes_to_CK^{lo}ER^{med}"                   
## [41] "Granulocytes_to_CK8-18^{+} ER^{hi}"                
## [42] "Granulocytes_to_MHC I^{hi}CD57^{+}"                
## [43] "HER2^{+}_to_CD57^{+}"                              
## [44] "HER2^{+}_to_CK8-18^{hi}ER^{lo}"                    
## [45] "HER2^{+}_to_MHC I^{hi}CD57^{+}"                    
## [46] "Ki67^{+}_to_Basal"                                 
## [47] "Ki67^{+}_to_CD57^{+}"                              
## [48] "Ki67^{+}_to_CK^{lo}ER^{med}"                       
## [49] "Ki67^{+}_to_CK8-18^{hi}ER^{lo}"                    
## [50] "Ki67^{+}_to_Ep CD57^{+}"                           
## [51] "Macrophages & granulocytes_to_CD57^{+}"            
## [52] "Macrophages & granulocytes_to_Ep CD57^{+}"         
## [53] "Macrophages_to_CD57^{+}"                           
## [54] "MHC I^{hi}CD57^{+}_to_CK8-18^{hi}ER^{lo}"          
## [55] "MHC I^{hi}CD57^{+}_to_Ep Ki67^{+}"                 
## [56] "MHC I^{hi}CD57^{+}_to_Fibroblasts FSP1^{+}"        
## [57] "MHC I^{hi}CD57^{+}_to_Granulocytes"                
## [58] "MHC I^{hi}CD57^{+}_to_HER2^{+}"                    
## [59] "MHC I^{hi}CD57^{+}_to_MHC^{hi}CD15^{+}"            
## [60] "MHC^{hi}CD15^{+}_to_B cells"                       
## [61] "MHC^{hi}CD15^{+}_to_CD57^{+}"                      
## [62] "MHC^{hi}CD15^{+}_to_CK^{+} CXCL12^{+}"             
## [63] "MHC^{hi}CD15^{+}_to_CK8-18^{+} ER^{hi}"            
## [64] "MHC^{hi}CD15^{+}_to_CK8-18^{hi}ER^{lo}"            
## [65] "MHC^{hi}CD15^{+}_to_Ep CD57^{+}"                   
## [66] "MHC^{hi}CD15^{+}_to_ER^{hi}CXCL12^{+}"             
## [67] "MHC^{hi}CD15^{+}_to_Fibroblasts FSP1^{+}"          
## [68] "MHC^{hi}CD15^{+}_to_MHC I^{hi}CD57^{+}"            
## [69] "T_{Reg} & T_{Ex}_to_CK8-18^{hi}ER^{lo}"

Investigate the cell counts of the combinations that could not be estimated

cells <- getCells()

# Barplot of cell type occurences in combinations
not_estimated_combinations_df <- tibble(combo = not_estimated_combinations) %>%
  separate(combo, into=c('phenotype_from','phenotype_to'),sep='_to_', remove=FALSE)

not_estimated_counts <- rbind(not_estimated_combinations_df %>% count(phenotype_from) %>% rename(phenotype = phenotype_from) %>% mutate(type = 'from'), 
                              not_estimated_combinations_df %>% count(phenotype_to)%>% rename(phenotype = phenotype_to) %>% mutate(type = 'to'))

not_estimated_counts_all <- merge(not_estimated_counts %>% expand(phenotype, type), not_estimated_counts, all.x=TRUE) %>%
                  mutate(n = coalesce(n, 0))
                              


ggplot(data=not_estimated_counts_all) + geom_bar(aes(x=reorder(phenotype, n), y=n, fill=type), stat='identity', position = position_dodge(), alpha = 0.75) + theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + ggtitle('Cell type occurences in combinations that could not be estimated') + ylab('n')

# Cell counts per sample
cell_counts <- cells %>% select(c(ImageNumber, meta_description)) %>%
                count(ImageNumber, meta_description) %>%
                rename(tnumber=ImageNumber) %>% mutate(type = ifelse(meta_description %in% not_estimated_combinations_df$phenotype_from &
                                                                       meta_description %in% not_estimated_combinations_df$phenotype_to, 'from & to', 
                                                                     ifelse(meta_description %in% not_estimated_combinations_df$phenotype_from, 'from',
                                                                      ifelse(meta_description %in% not_estimated_combinations_df$phenotype_to, 'to','not'))))

group_ordered <- with(cell_counts, reorder(meta_description,n,median))

cell_counts_ordered <- cell_counts
cell_counts_ordered$meta_description <- factor(cell_counts_ordered$meta_description,levels = levels(group_ordered))


# a <- ifelse(unique(cell_counts$meta_description) %in% unique(not_estimated_combinations_df$phenotype_from), "red", "black")

ggplot(cell_counts_ordered, aes(x=meta_description,n,y=n, fill=type)) + 
  geom_boxplot(outlier.shape = NA) + theme_bw() +
  scale_fill_manual(values=c("darkgrey", 'red', 'white',"lightgrey")) +
  scale_y_continuous(limits = quantile(cell_counts_ordered$n, c(0.1, 0.9))) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + ggtitle('Counts per cell type')
## Warning: Removed 1524 rows containing non-finite values (`stat_boxplot()`).

###
cell_occurences <- merge(cell_counts %>% expand(tnumber, meta_description), cell_counts, by = c('tnumber','meta_description'), all.x = TRUE) %>%
                  mutate(n = coalesce(n, 0)) %>%
                  mutate(category = ifelse(n==0,'0', ifelse(n<10, '0<10', ifelse(n<50,'10<50','>50')))) %>% group_by(meta_description) %>%
                  count(category)


lvls <- cell_occurences %>% filter(category == '0') %>% arrange(desc(n)) %>% pull(meta_description)
a <- ifelse(lvls %in% unique(not_estimated_combinations_df$phenotype_from), "red", "black")

ggplot(cell_occurences, aes(x=factor(meta_description, level = lvls), y=n, fill=factor(category, levels=c('>50','10<50','0<10','0')))) + 
    geom_bar(position="stack", stat="identity") + theme_bw()  +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, colour=a)) + ggtitle('Cell type occurences in samples') + ylab('number of samples') +   xlab('cell type')+ guides(fill=guide_legend(title="occurences"))
## Warning: Vectorized input to `element_text()` is not officially supported.
## ℹ Results may be unexpected or may change in future versions of ggplot2.

Initial parameters could be fitted so we can have a look at these estimations.

ggplot(data = initial_parameters, aes(x=shape, y=scale)) +
  geom_point(alpha=0.5, size=0.5) +
  geom_point(data = initial_parameters %>% filter(phenotype_from == 'CD57^{+}'), aes(x=shape, y=scale,color=combo)) +
  theme(legend.position='right') +
  scale_y_log10() + theme_bw() +
  ylim(c(10,300)) + xlim(0,6) +
  xlab('Shape') + ylab('Scale')
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Warning: Removed 1360 rows containing missing values (`geom_point()`).
## Warning: Removed 275 rows containing missing values (`geom_point()`).

ggplot(data = initial_parameters, aes(x=shape, y=scale)) +
  geom_point(alpha=0.5, size=0.5) +
  geom_point(data = initial_parameters %>% filter(phenotype_to == 'CD57^{+}'), aes(x=shape, y=scale,color=combo)) +
  theme(legend.position='right') +
  scale_y_log10() + theme_bw() +
  ylim(c(10,300)) + xlim(0,6) +
  xlab('Shape') + ylab('Scale')
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Warning: Removed 1360 rows containing missing values (`geom_point()`).
## Warning: Removed 94 rows containing missing values (`geom_point()`).

plot_initial_parameters = function(t,c){
  shape_param = filter(initial_parameters, tnumber==t & combo==c)['shape']
  scale_param = filter(initial_parameters, tnumber==t & combo==c)['scale']
  print(ggplot(all_distances_data[all_distances_data$tnumber == t & all_distances_data$phenotype_combo == c ,]) +
           geom_bar(aes(x=distance_window, y=N.per.mm2.scaled), stat="identity") + theme_bw() +
           stat_function(fun = dweibull, args = list(shape = shape_param,
                                                     scale = scale_param))
           + labs(title = paste("tnumber: ", t, "combi: ", c),
              subtitle = paste('parameters: ', shape_param, ' ', scale_param)) + xlim(0,300))
}

combination = not_estimated_combinations[1]
plts_notexistingparams <- lapply(sample(unique(all_distances_data %>% filter(phenotype_combo == combination) %>% pull(tnumber)), 3, replace=FALSE), plot_initial_parameters, c=not_estimated_combinations[1])
## Warning: Computation failed in `stat_function()`
## Caused by error in `fun()`:
## ! Non-numeric argument to mathematical function

## Warning: Computation failed in `stat_function()`
## Caused by error in `fun()`:
## ! Non-numeric argument to mathematical function

## Warning: Computation failed in `stat_function()`
## Caused by error in `fun()`:
## ! Non-numeric argument to mathematical function

# gridExtra::grid.arrange(grobs = plts_notexistingparams, heights=c(10,10,10))

Parameter plots

Do the parameters look like a C-curve?

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Warning: Removed 107 rows containing non-finite values (`stat_density2d()`).
## Warning: The following aesthetics were dropped during statistical transformation: label
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?
## Warning: Removed 107 rows containing missing values (`geom_point()`).

## Warning in geom_point(data = means %>% mutate(type = "means"), alpha = 0.5, :
## Ignoring unknown aesthetics: label
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Warning: Removed 107 rows containing missing values (`geom_point()`).

## Rows: 1152 Columns: 11
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (6): sample, phenotype_combo, translational_response, Biopsy, PA_Respons...
## dbl (5): a, b, A, B, ptID
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Warning: Removed 107 rows containing missing values (`geom_point()`).
## Warning: Removed 1 rows containing missing values (`geom_point()`).

Make histograms of curve extremes

# Pick extremes
# Divide the curve in areas and pick random from these areas
bottom_right <- sample_n(all_parameters %>% filter(shape > 5) %>% filter(scale < 20),1) %>% mutate(type= "low scale, high shape")
bottom_left <- sample_n(all_parameters %>% filter(shape < 1.5) %>% filter(scale < 20),1)%>% mutate(type= "low scale, low shape")
top_right <- sample_n(all_parameters %>% filter(shape > 3) %>% filter(scale > 300),1)%>% mutate(type= "high scale, high shape")
top_left <- sample_n(all_parameters %>% filter(shape < 3) %>% filter(scale > 300),1)%>% mutate(type= "high scale, low shape")
belly <- sample_n(all_parameters %>% filter(shape < 3) %>% filter(shape > 1) %>% filter(scale > 50) %>% filter(scale <60),1)%>% mutate(type= "belly")

sample_points <- rbind(bottom_right, bottom_left, top_right, top_left, belly) %>% mutate(color = c('red', 'blue', 'yellow', 'green', 'purple'))

ggplot() +
  geom_point(data = all_parameters, aes(x=shape, y=scale), alpha=0.5, size=0.5) +
  geom_point(data = sample_points, aes(x=shape, y=scale), alpha=1, size=3, color = sample_points$color) +
  theme_bw() +
  ylim(10,300) + xlim(0,6) +
  scale_y_log10() +
  xlab('Shape') + ylab('Scale') + ggtitle("parameters with densities")
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Warning: Removed 107 rows containing missing values (`geom_point()`).

ggplot() + geom_bar(data  = all_distances_data[all_distances_data$tnumber == bottom_left$tnumber & 
                                        all_distances_data$phenotype_combo == bottom_left$phenotype_combo,], 
           aes(x=distance_window, y=N.per.mm2.scaled, color=bottom_left$type), stat="identity") +
         stat_function(fun = dweibull, args = list(shape = bottom_left$shape, scale = bottom_left$scale)) + 
  # geom_bar(data  = all_distances_data[all_distances_data$tnumber == bottom_right$tnumber & 
  #                                     all_distances_data$phenotype_combo == bottom_right$phenotype_combo,], 
  #        aes(x=distance_window, y=N.per.mm2.scaled), stat="identity") +
  #        stat_function(fun = dweibull, args = list(shape = bottom_right$shape, scale = bottom_right$scale), colour='green') + 
  # geom_bar(data  = all_distances_data[all_distances_data$tnumber == top_left$tnumber & 
  #                                       all_distances_data$phenotype_combo == top_left$phenotype_combo,], 
  #          aes(x=distance_window, y=N.per.mm2.scaled), stat="identity") +
  #          stat_function(fun = dweibull, args = list(shape = top_left$shape, scale = top_left$scale), colour='purple') + 
  # geom_bar(data  = all_distances_data[all_distances_data$tnumber == top_right$tnumber & 
  #                                       all_distances_data$phenotype_combo == top_right$phenotype_combo,], 
  #          aes(x=distance_window, y=N.per.mm2.scaled), stat="identity") +
  #          stat_function(fun = dweibull, args = list(shape = top_right$shape, scale = top_right$scale), colour='red') + 
  xlim(0,300)

Now look at specific ‘X to celltype’ parameters.

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Warning: Removed 3 rows containing missing values (`geom_point()`).
## Saving 7 x 5 in image
## Warning: Removed 3 rows containing missing values (`geom_point()`).
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

## Saving 7 x 5 in image
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

## Saving 7 x 5 in image
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

## Saving 7 x 5 in image
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

## Saving 7 x 5 in image
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

## Saving 7 x 5 in image

Heatmap

What patterns in the parameters do we see?

#Generate a correlation matrix

parameter_matrix <- tibble(tnumber = unique(all_parameters$tnumber))

for (c in unique(all_parameters$phenotype_combo)){
  parameter_matrix <- left_join(x = parameter_matrix, y= all_parameters %>% 
                              select(c(tnumber, phenotype_combo, shape)) %>% 
                              filter(phenotype_combo == c),
                            by='tnumber') %>%
                      select(-c(phenotype_combo))
  names(parameter_matrix)[names(parameter_matrix) == 'shape'] = c
}

parameter_matrix <- parameter_matrix %>% select(-c(tnumber))

correlation_matrix <- Hmisc::rcorr(as.matrix(parameter_matrix), type="pearson")
## Warning in sqrt(npair - 2): NaNs produced
## Warning in pt(abs(h) * sqrt(npair - 2)/sqrt(pmax(1 - h * h, 0)), npair - : NaNs
## produced
flattenCorrMatrix <- function(cormat, pmat, nmat) {
  ut <- upper.tri(cormat)
  data.frame(
    row = rownames(cormat)[row(cormat)[ut]],
    column = rownames(cormat)[col(cormat)[ut]],
    cor  =(cormat)[ut],
    p = pmat[ut],
    n = nmat[ut]
    )
}

flattened_correlation_matrix <- flattenCorrMatrix(correlation_matrix$r, correlation_matrix$P, correlation_matrix$n) %>% drop_na()  %>%
  filter(n > 200) %>%
  mutate(sig_p = ifelse(p < .05, T, F), p_if_sig = ifelse(p <.05, p, NA), r_if_sig = ifelse(p <.05, cor, NA))

flattened_correlation_matrix %>% 
 ggplot(aes(row, column, fill=cor, label=round(r_if_sig,2))) +
 geom_tile() + 
 labs(x = NULL, y = NULL, fill = "Pearson's\nCorrelation")+ 
  theme(axis.text = element_text(size = 1))  + 
 scale_fill_gradient2(mid="#FBFEF9",low="#0C6291",high="#A63446", limits=c(-1,1)) +
 scale_x_discrete(expand=c(0,0)) + 
 scale_y_discrete(expand=c(0,0))

significant_correlations <- flattened_correlation_matrix %>% 
  filter(sig_p == TRUE) %>%
  filter(p_if_sig > 0) %>%
  arrange(p_if_sig) 

significant_correlations %>% head()
##                                         row
## 1                Fibroblasts_to_Fibroblasts
## 2            CK^{med}ER^{lo}_to_Macrophages
## 3             Myofibroblasts_to_Fibroblasts
## 4          Myofibroblasts_to_CK^{lo}ER^{lo}
## 5                Fibroblasts_to_Macrophages
## 6 CK^{med}ER^{lo}_to_CK8-18^{hi}CXCL12^{hi}
##                                     column       cor            p   n sig_p
## 1           CK^{med}ER^{lo}_to_Fibroblasts 0.3500501 2.220446e-16 520  TRUE
## 2         CK^{+} CXCL12^{+}_to_Macrophages 0.5345669 2.220446e-16 204  TRUE
## 3           CK^{med}ER^{lo}_to_Fibroblasts 0.4193625 4.440892e-16 345  TRUE
## 4           Fibroblasts_to_CK^{med}ER^{lo} 0.3839297 4.440892e-16 421  TRUE
## 5   Myofibroblasts PDPN^{+}_to_Macrophages 0.4919914 4.440892e-16 241  TRUE
## 6 CK^{lo}ER^{lo}_to_CK8-18^{hi}CXCL12^{hi} 0.5100846 4.440892e-16 220  TRUE
##       p_if_sig  r_if_sig
## 1 2.220446e-16 0.3500501
## 2 2.220446e-16 0.5345669
## 3 4.440892e-16 0.4193625
## 4 4.440892e-16 0.3839297
## 5 4.440892e-16 0.4919914
## 6 4.440892e-16 0.5100846
for (c in 1:5){
point <- all_parameters %>%
    filter(phenotype_combo == significant_correlations[c,"row"] | phenotype_combo == significant_correlations[c,"column"]) %>%
    ggplot(aes(x=shape, y=scale, color=phenotype_combo)) +
    geom_point(alpha=0.4, size=3) +
    ylim(10,300) + xlim(0,6) +
    theme_bw() + scale_y_log10() +
    xlab('Shape') + ylab('Scale') +
    theme(legend.position='right') + ggtitle(paste('correlated with p-value: ', signif(significant_correlations[c,"p_if_sig"],digits=3))) + 
  guides(colour = guide_legend(override.aes = list(size=10)))
  
  print(point)
}
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

## Warning: Removed 1 rows containing missing values (`geom_point()`).
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.